SCIENTIFIC 

REPORTS 




open The Boson peak in supercooled water 



SUBJECT AREAS: 

PHASE TRANSITIONS 
AND CRITICAL 
PHENOMENA 

FLUIDS 

THERMODYNAMICS 

STRUCTURE OF SOLIDS AND 
LIQUIDS 



Received 
26 November 2012 



Accepted 
22 May 2013 

Published 
17 June 2013 



Correspondence and 
requests for materials 
should be addressed to 
P.K. (pradeep.kumar® 
rockefeller.edu] or 
K.T.W. (wikfeldt@hi.is) 



Pradeep Kumar', K. Thor Wikfeldt 23 , Daniel Schlesinger 4 , Lars G. M. Pettersson 4 & H. Eugene Stanley 5 

'Center for Studies in Physics and Biology, The Rockefeller University, 1 230 York Avenue, New York, NY 1 0021 USA, 2 Science 
Institute and Faculty of Science, VR-III, University of Iceland, 1 07 Reykjavik, ICELAND, 3 Nordifa, Royal Institute of Technology and 
Stockholm University, Roslagstullsbacken 23, S-l 0691 Stockholm, SWEDEN, 4 Department of Physics, AlbaNova, Stockholm 
University, S-l 0691 Stockholm, SWEDEN, 5 Center for Polymer Studies and Department of Physics, Boston University, Boston, MA 
02215 USA. 



We perform extensive molecular dynamics simulations of the TIP4P/2005 model of water to investigate the 
origin of the Boson peak reported in experiments on supercooled water in nanoconfined pores, and in 
hydration water around proteins. We find that the onset of the Boson peak in supercooled bulk water 
coincides with the crossover to a predominantly low-density-like liquid below the Widom line T w . The 
frequency and onset temperature of the Boson peak in our simulations of bulk water agree well with the 
results from experiments on nanoconfined water. Our results suggest that the Boson peak in water is not an 
exclusive effect of confinement. We further find that, similar to other glass-forming liquids, the vibrational 
modes corresponding to the Boson peak are spatially extended and are related to transverse phonons found 
in the parent crystal, here ice Ih. 

One of the characteristic features of many glasses and amorphous materials is the onset 1 of low-frequency 
collective modes (Boson peak) in the energy range 2-10 meV at low T, where the vibrational density of 
states (VDOS) g(oj) shows an excess over g(uj) oc co 2 predicted by the Debye model. Disordered 
materials are further known to exhibit many anomalous behaviors compared to their crystalline counterparts, 
such as the temperature dependence of thermal conductivity 23 and specific heat 4,5 at low temperatures. Many 
scenarios 6 ' 7 have been suggested to explain the physical mechanisms behind the Boson peak and related anom- 
alies, but a comprehensive understanding has proved elusive. 

Recent neutron scattering experiments on water confined in nanopores indicate the presence of a Boson peak 8,9 
around 5-6 meV (40 - 49 cm" 1 ) emerging below 230 K in the incoherent dynamic structure factor. These results 
were tentatively interpreted as arising from a gradual change in the local structure of confined liquid water when 
crossing the Widom line temperature TV 10 . Earlier, neutron scattering has also been applied to protein hydration 
water 11 and a Boson peak was found around 30 cm" 1 . T w corresponds to the loci of maxima of thermodynamic 
response functions in the one-phase region beyond the liquid-liquid critical point (LLCP) proposed to exist in 
supercooled liquid water 12 . A Widom line in the liquid-gas supercritical region in argon has recently been 
studied 13,14 and found to be directly related to a dynamical crossover between liquid-like and gas-like properties, 
but the existence of a dynamical crossover in supercooled water is subject to some controversy 15 23 . Since the 
melting temperature is strongly depressed in nanoconfinement and in protein hydration water, deeply super- 
cooled confined water has been used experimentally to infer the behavior of bulk water, which is more challenging 
to supercool. However, the similarity to bulk water has been called into question 24,25 , and it is thus important to 
further investigate the relationship between water in these different forms. 

Both experimental 9,26 and simulation 27 studies suggest that the density relaxation of confined and hydration 
water at and slightly below the Widom temperature is of the order of a few tens of nanoseconds, implying that 
liquid water is still in metastable equilibrium over the experimental time scales involved 28 . The experimental 
observation of a Boson peak below T w in confined water thus suggests that the low-density-like liquid shares 
vibrational properties with the glassy state, as has been observed previously for other systems such as B 2 0 3 29,30 . On 
the other hand, the dynamic structure factor of crystalline ice exhibits a peak at a slightly higher frequency around 
50 cm 4 3133 as do the Raman spectra of ice Ih and proton-ordered ice XI 34 ; in the latter case the peak becomes 
extremely sharp. These results indicate a connection between the dynamics of supercooled liquid and crystalline 
vibrational dynamics. Indeed, it has recently been suggested that Boson peaks observed in glasses are related to the 
vibrational dynamics of the parent crystal 35 . 

To shed light on the question of whether bulk supercooled liquid water displays a Boson peak, we study the low- 
frequency dynamics of the TIP4P/2005 model of water which accurately reproduces a range of water properties 36 , 
including its anomalies 37 . This model has been found to exhibit a LLCP in the vicinity of P c = 1350 bar and T c = 
193 K 38 . The associated Widom line has been shown to be accompanied by a structural crossover from a 
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predominantly high-density liquid (HDL) structure to a predomi- 
nantly low-density liquid (LDL) structure 39,40 occurring at a temper- 
ature T w around 230 K at 1 bar. We find that this water model 
indeed displays a Boson peak in the bulk supercooled regime and 
that its onset coincides with the Widom temperature. Our analysis 
further shows that it derives from transverse acoustic modes in the 
parent crystal, ice Ih, in agreement with the recently proposed pic- 
ture 35 . We further verify our results using another model of water, 
TIP5P, and find also in this case an emergence of a Boson peak below 
the Widom line. 

Results 

Incoherent dynamic structure factors. A quantity that is readily 
accessible in inelastic neutron scattering experiments is the 
incoherent dynamic structure factor (DSF) Ss(k,oj) which probes 
single-particle dynamics. To compare our simulation results with 
experiment we first calculate the self- intermediate scattering 
function 



v[T,(t)-r,{0)] 



(1) 



where rj(f) are the positions of oxygen atoms, N is the number of 
molecules and angular brackets denote an ensemble average and 



averaging over different k with the same modulus. In simulations 
the wave vector k is defined as ^ (n x , n y , n z ) for integers (n x , n y , n z ) 
and system size L. We perform the frequency decomposition of 
Fs(k,t) to obtain the incoherent DSF 



S s (k, w) 



F s (k, t)e WJ 'dt. 



(2) 



In Figs. 1(a) and 1(b), we show Fs(k,t) and Ss(k,co) for supercooled 
liquid water simulations with N = 512 molecules at atmospheric 
pressure and different temperatures for k = 2 A" 1 , i.e., the position 
of the first peak in the structure factor S(k). For T < T w , a minimum 
appears in F s (k,t) around 0.3 ps followed by oscillations up to around 
10 ps. In the incoherent DSF these changes with temperature corre- 
spond to the emergence of sharp peaks at 20 and 25 cm -1 along with a 
broad peak centered around 37 cm -1 . For comparison, we also show 
F s (k,t) and S s (k,co) from hexagonal ice simulations at 100 K for two 
different system sizes. The connection between liquid and crystalline 
low- frequency dynamics will be further discussed below but we note 
here that the system size obviously affects the region co < 40 cm" 1 . 

Since low-frequency dynamics in supercooled liquids and glasses 
may be affected by finite-size effects 41-43 , we show in Figs. 1(c) and 
1(d) a comparison with a much larger system, N = 45, 000 at 210 K. 
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Figure 1 | (a) F s (k,t) and (b) S s (k,m) for a range of temperatures 210-260 K at fixed k = 2.0 A -1 and system size JV = 512. S s (k,oj) at k = 2.0 A -1 for ice at 
100 K is also shown for two different system sizes, N = 432 and JV = 3456. (c)-(d) show a comparison with F s (k,t) and Ss(k,co) for a large N = 45,000 
simulation at 210 K at k ~ 1.6 A -1 . The broad peak around co = 37 cm -1 is evidently independent of system size while the sharper peaks at lower 
frequency are not present in the large simulation. 
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In F s (k,t) the oscillations are shifted to longer times while the min- 
imum at 0.3 - 0.4 ps and subsequent peak around 0.9 ps are system- 
size independent. The sharp peaks in Ss(k,oS) vanish in the large 
simulation but the broad peak around 37 cm" 1 persists; we therefore 
label this peak of low-frequency excitations as the Boson peak. 
Comparing our simulation results with experimental neutron data 
on protein hydration water 11 and on water in confinement 8,9 reveals 
rather good agreement. Both the experimental energy position, 
30 cm" 1 for hydration water 11 and 45 cm" 1 for confined water 8,9 , 
and the temperature onset, T = 225 K 8,9 , are surprisingly well repro- 
duced by the TIP4P/2005 model considering the approximate nature 
of classical force fields. Since we observe well-defined low-frequency 
modes in the incoherent DSF in the simulations of bulk water, we 
conclude that the Boson peak in supercooled water is not a con- 
sequence of confinement and that it would likely be detected also 
in experiments on bulk water, if sufficient supercooling could be 
achieved. 

A connection to the Widom line in TIP4P/2005 is clearly present 
since a qualitative change of both Fs(k,t) and Ss(k,oj) takes place 
between 230 and 240 K. Our results thus imply that a Boson peak 
may also appear in other tetrahedral 44 liquids, such as silicon and 
silica, for which simulations have indicated the existence of a liquid- 
liquid phase transition 45,46 . To confirm the connection between the 
Boson peak and the Widom line we have performed additional simu- 
lations of a different water potential, the TIP5P model 47 , which has 
been shown to exhibit a Widom line around 250 K 27 , 20 K above that 
for TIP4P/2005. As shown in the Supplementary Information, a 
Boson peak emerges also in this model and, indeed, it coincides with 
T w at 250 K. 

Vibrational density of states. The Boson peak is commonly 
discussed in terms of the VDOS, g(to). A distinct peak is often not 
seen in the VDOS itself but appears only in the reduced VDOS after 
normalizing by the squared frequency, g(co)/co 2 , which reveals the 
excess over the Debye model prediction g(w) oc co 2 . We calculate 
g(oj) as the Fourier transform of the oxygen velocity autocorrelation 
function, C v (f): 



where 



c v« = (^I>(fH-(°) 



(4) 



C v {i)e ,mt dt, 



(3) 



The sum includes all N oxygen atoms in the system, V; are the oxygen 
atom velocities and (...) denotes the ensemble average. 

In Fig. 2(a) we show graphs of g(co) for various temperatures. 
Below T w ,g{co) shows an onset of the same two sharp low- frequency 
peaks as observed in Ss(k,co) in Fig. 1. By simulating a range of 
different system sizes we can establish the system-size dependence 
of g(oj). The agreement for frequencies co > 40 cm" 1 is very good for 
different system sizes L, but the sharp low-frequency peaks shift to 
lower frequencies as the system size is increased. Indeed, both peaks 
extrapolate linearly to zero as 1/L, as seen in the inset of Fig. 2(a), 
suggesting that they disappear in the limit of infinite system size. 

As discussed in relation to Fig. 1, there is a low-frequency peak in 
Ss(k,w) for hexagonal ice at somewhat higher frequency compared to 
the supercooled liquid simulations, suggesting a link between the 
liquid and crystalline low-frequency modes. We investigate this by 
performing thermally non-equilibrated simulations at even lower 
temperatures where the TIP4P/2005 model vitrifies to a low-density 
amorphous (LDA)-like solid. The non-equilibrium simulations are 
performed by annealing the equilibrated metastable 210 K simu- 
lation with a cooling rate of 2 - 10 10 K/s to reach target temperatures 
of 150 and 100 K. Figure 2(b) shows the reduced VDOS, lg(co) - 
g(0)]/co 2 for simulations at 100, 150 and 210 K compared to the 
hexagonal ice simulation at 100 K. We subtract an extrapolated value 
ofg(0) iromg(cu) to eliminate uncertainties due to the finite length of 
the trajectories when evaluating the Fourier transform in Eq. (3). A 
clear shift to higher frequencies of the Boson peak is observed as the 
supercooled liquid simulations are cooled into the LDA glass region, 
and the reduced VDOS at 100 K resembles the crystalline ice 
counterpart. 

In the inset of Fig. 2(b), we show the reduced VDOS gtaA.co)lco 2 
obtained from the normal modes of the liquid calculated from 
quenched configurations (i.e., energy minimized inherent struc- 
tures) at T = 210 K. We find that the Boson peak in gNM^AIoi 2 is 
blue-shifted to around 50 cm" 1 , close to the peak for hexagonal ice, 
compared to the velocity autocorrelation function VDOS which 
peaks at 37 cm" 1 . This suggests that a key difference in low- frequency 




Figure 2 | (a) VDOS of liquid water at supercooled temperatures and of hexagonal ice and (inset) inverse box-length dependence of the two lowest sharp 
peaks of the VDOS, extrapolating to zero frequency in the limit of infinite system size. In (a) the lowest-fc transverse current spectrum is shown in the 
bottom part, illustrating that the sharp low-frequency peaks are low-fc transverse modes, (b) Reduced vibrational density of states at low Tcalculated as 
[g(co) — g(0)]/co 2 with an extrapolated g(0) to eliminate uncertainties related to the finite simulation time. Upon rapid cooling into a non-equilibrated 
LDA ice at 150 K and 100 K, the Boson peak is seen to shift to higher frequencies, approaching that of hexagonal ice. The inset in (b) shows the reduced 
VDOS calculated from the normal modes of inherent structures quenched from equilibrated T = 210 K configurations. The normal mode gjvM^)/™ 2 
shows a Boson peak which is blue-shifted to about the same frequency as the crystal, suggesting that as the liquid structure is made harmonic the Boson 
peak frequency shifts to higher values saturating at about 50 cm" 1 . 
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vibrational properties between the LDL-like liquid below the Widom 
line and crystalline ice lies in the more anharmonic dynamics of the 
liquid phase. 

Transverse and longitudinal correlation functions. Having 
established the existence of a Boson peak in supercooled water 
below the Widom line and its connection to low-frequency 
dynamics present in the parent crystal, ice Ih, we now turn to the 
study of transverse and longitudinal current correlations to clarify 
the nature of these low-frequency modes. We calculate longitudinal 
and transverse currents as 



/i(k, r)=^k r v i (r)exp[-tkT j ( 



7r(k, r)=^k x -v j (r)exp[-tkT j (t) 



(5) 



(6) 



where k|| and ki are unit vectors respectively parallel and perpen- 
dicular to k, and r,(r) and v,(f) denote the oxygen atoms' position and 
velocity, respectively. The frequency decomposition of the longitu- 
dinal and transverse current autocorrelation functions is 



Q(/C, Co) 



(L(k, t)U-k, 0))e io *dt 



where a = L or T. 

In the bottom part of Fig. 2(a) we show superimposed on the 
VDOS the transverse current correlation function (TCCF) Cj{k,io) 
at different temperatures for the lowest wave number k accessible in 
the N = 512 simulation boxes, i.e., k = 2rc(l, 0, 0)/L and permuta- 
tions thereof. We see that the first sharp size-dependent low- 
frequency peak in the VDOS, which develops below T w , coincides 
exactly with the lowest-fc TCCF, and the second finite-size peak 
around 25 cm" 1 coincides with the second-lowest-A: TCCF (not 
shown). Returning to Fig. 1(d), the low-frequency side of the 
Boson peak is smoother in the large N = 45, 000 simulations and 
the sharp peaks related to the lowest-fc transverse currents are instead 
seen at frequencies below 10 cm" 1 . We thus conclude that the sharp 
system size dependent low-frequency peaks around the Boson peak 
frequency and below are predominantly transverse excitations, con- 
sistent with previous findings for amorphous silica 43 . 

The dispersion relations for transverse and longitudinal current 
spectra, Cj{k,co) and Ci(k,w) at T = 210 K, are shown in Fig. 3 (the 
current spectra are shown in the Supplementary Material). We 
obtain dispersion relations by fitting damped harmonic oscillator 
(DHO) 48 lines to both longitudinal and transverse spectra. 

At small k below 0.5 A -1 , the longitudinal current spectrum, 
Q.(fc,co), shows only one acoustic dispersing excitation. For k > 0.5 
A -1 , Ci(k,co) shows the existence of three excitations and can be fit 
with three DHO lines. Besides the dispersing excitation at intermedi- 
ate frequency, two other weakly dispersing excitations appear — one 
at low frequency around 50 cm" 1 and the other at high frequency 
around 260 cm" 1 . The intensity of these excitations in C L (k,co) 
increases upon further increase of k. Transverse current spectra 
Cj{k,co) exhibit an acoustic dispersing excitation for k < 0.5 A" 1 . 
For k > 0.5 A" 1 , C r (fc,co) develops a peak at co ~ 260 cm" 1 , the same 
excitation as in the longitudinal current spectra. Transverse and 
longitudinal modes are thus strongly mixed above 0.5 A" 1 as has been 
found previously also for the SPC/E water model 49 . Moreover, the 
mixing happens at both high frequencies (co ~ 260 cm" 1 ) and low 
frequencies (co < 60 cm" 1 ) as evident from transverse excitations 
appearing in longitudinal spectra. We note that the band at 
260 cm" 1 in liquid water is associated with four-coordinated water 
molecules since low-frequency Raman spectra of water down to - 
20°C showed its intensity to increase with decreasing temperature 50 
and hexagonal ice also displays a strong band near this frequency 51 . 




Figure 3 | Longitudinal (filled circles) and transverse (filled squares) 
dispersion relations calculated from the peak positions of co L ( k) in C L { fc,co) 
and cojik) in C T (fc,to) for N = 2048 and T = 210 K. For k < 0.5 A" 1 , 
longitudinal spectra exhibit one dispersive branch (LA: black filled circles), 
while for k > 0.5 A" 1 , the longitudinal spectra exhibit three excitations - 
low and high frequency weakly dispersive excitations (shown in orange and 
red filled circles). For k < 0.5 A" 1 , transverse current spectra exhibit only 
one dispersive excitation (TA: blue filled squares) while for k > 0.5 A -1 , it 
exhibits both the dispersive branch as well as a weakly dispersive excitation 
(green filled squares) at co ~ 260 cm" 1 . For information on extraction of 
the dispersion relation see Fig. S2. 

Hence, the emergence of the 260 cm" 1 band in both C L (k,co) and 
Cj{k,w) suggests that liquid water at this temperature exhibits net- 
works of four-coordinated molecules over length scales as large as 
27t/0.5 A" 1 ~ 13 A. This observation is also consistent with recent 
studies where it is shown that the sizes of clusters of highly tetrahed- 
ral molecules increases below the Widom temperature 3,44 . 

The emergence of an additional, high-frequency excitation in 
C L (k,co) and C r (fc,co) around 260 cm" 1 for k > 0.5 A" 1 suggests a 
low-frequency liquid-like and a high-frequency solid-like behavior 
of the longitudinal and transverse spectra at this length scale and a 
concomitant pile-up of spectral intensity takes place in the Boson 
peak regime. 

To get further insights into the longitudinal or transverse char- 
acter of the Boson peak, we compare cumulatively integrated spectra 
over k for different frequencies for both the transverse and longit- 
udinal parts, C L (k, co) and Cj(k, co), defined as 



C*(ha>) = 



C a (k! , co)dk', 



(8) 



where a = L or T and k„ 



2n/L is the smallest wave number 



accessible in our system of box size L. C*(k, co) thus describes the 
total contribution of longitudinal and transverse modes with differ- 
ent wave numbers up to k for a given frequency co. In Fig. 4(a)-(b) we 
show normal and cumulatively integrated current spectra for several 
frequencies co B in the Boson peak region around 37 cm" 1 . It can be 
clearly seen that transverse modes are dominant in the Boson peak 
frequency regime for k > 0.5 A" 1 . 

Localization analysis. A number of studies on amorphous materials 
have found that the modes in the Boson peak frequency range are 
localized or quasi-localized 52,53 . In order to investigate whether this 
holds also for TIP4P/2005 water below the Widom line T w , we 
calculate the degree of localization by performing a normal mode 
analysis of quenched (or inherent) structures, obtained by energy- 
minimizing snapshots from simulation trajectories. A measure of the 
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Figure 4 | (a) Current correlation functions C L {k,to B ) and C-j{k,co B ) as function of k. (b) Integrated current correlation functions C* L (k, co B ) and 
Cy(fc, (Ob) as function of k (see Eq. 8). Several frequencies to B in the Boson peak frequency regime are shown. 



degree of localization of a vibrational mode is the frequency- 
dependent participation ratio 54,55 p fl for mode [i 

i -i 



(9) 



where vt is the contribution of all degrees of freedom of molecule i to 
the normal mode ji. The participation ratio is unity when all 
molecules contribute equally to the normal mode in consideration, 
while p = l/Nii only one molecule contributes to the total energy of 
the mode. Hence, for an extended mode^, is close to unity and does 
not depend on the system size while for a localized mode it is small 
and scales with system size as l/N. One way to determine the extent 
of localization is thus to inspect the system-size dependence of the 
participation ratio. We find (see Fig. 5(a)) that the participation ratio 
is quite large, around 0.6, for the modes with frequency below 50 cm" 
the region of the Boson peak. Apart from the sharp finite-size 
peaks, the participation ratio for the modes in the region of the 
Boson peak show only a weak system size dependence, suggesting 
that the modes giving rise to the Boson peak are not localized but 
extended. 

We next introduce a function A max (r) defined as the maximum 
displacement of molecules at a distance r from the molecule with the 
largest displacement in the normal mode. For a localized mode a 



rapid decay of A max (r) should be seen. In Fig. 5(b) we show average 
A max (r) for the normal modes in two different frequency regimes - 
for the Boson peak regime (co B = 45 ± 2.5 cm -1 , note that the 
frequency of the Boson peak for the quenched configuration is 
shifted to higher frequency compared to the Boson peak frequency 
co B ~ 37 cm -1 of the liquid at T = 210 K, see inset of fig. 2(b)), and in 
comparison we also show average A max (r) for modes in the range of 
co = 400 ± 5 cm 4 , which is a range of localized vibrations. While for 
co = 400 ± 5 cm" 1 , A max (r) decays rapidly to zero, for the Boson peak 
region it does not, suggesting again an extended nature of the Boson 
peak modes. 

Discussion 

The low-frequency vibrations of a classical potential model of water, 
TIP4P/2005, are investigated in the supercooled temperature regime 
to clarify the origin of the Boson peak reported from inelastic neutron 
scattering experiments below around 225 K in nano-confmed 8,9 and 
protein-hydration water 11 . We find that sharp low-frequency peaks 
emerge in the incoherent dynamic structure factor and the reduced 
density of states of the simulated liquid water as the system is cooled 
below the Widom line, but a system-size investigation reveals that in 
the limit of an infinitely large simulation box these peaks extrapolate 
to zero frequency. The sharp finite-size peaks are seen to coincide 
exactly with the low wave-number acoustic branch of the transverse 




Figure 5 | (a) Frequency dependent participation ratio for two system sizes, N= 512 and N = 1024. (b) Spatial dependence of the amplitude of normal 
modes A max (r) in the Boson peak regime. For a comparison, A max (r) for co = 400 ± 5 cm" 1 where the modes are localized is also shown. The value of 
A max (0) is normalized to 1 for the range of frequencies shown. 
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current correlation functions, reflecting a strong contribution of 
transverse modes in this frequency region. However, we find a broad 
Boson peak centered around 37 cnr 1 which is unaffected by system 
size, and for which the frequency region and temperature of onset in 
the incoherent DSF agree well with neutron experiments on confined 
water. Due to its lower melting temperature, water in confinement 
has been used experimentally to infer the behavior of bulk water 
below the bulk homogeneous nucleation temperature. The validity 
of this comparison has been questioned, but the good agreement 
observed here in the low-frequency vibrational dynamics lends sup- 
port to the view of confinement as useful in the study of supercooled 
bulk water, at least for low-frequency vibrational properties. 

The frequency of the Boson peak in supercooled TIP4P/2005 
water as observed in the reduced VDOS changes as the simulation 
is annealed into the LDA glass region and approaches cu = 45 cm" 1 . 
This is also the frequency at which hexagonal ice simulations display 
a peak in vibrational spectra deriving from the transverse acoustic 
branch, as has been observed experimentally 31 33 . Thus, upon low- 
ering the temperature below the Widom line, the low-frequency 
dynamics of the system progressively changes from LDL-like to 
LDA glass and to the dynamics found in hexagonal ice. A similar 
shift to higher frequencies is observed in normal-mode spectra of 
inherent structures quenched from liquid at temperatures below the 
Widom line, indicating that the lower frequency of the Boson peak in 
the liquid below the Widom line, compared to the transverse acoustic 
peak in the ice, is a result of the more anharmonic nature of the 
vibrational modes in the liquid. 

Recent work by Chumakov et aZ. 35 on glasses suggests that there is 
no excess in the actual number of states at the Boson peak and hence 
no additional modes compared to the crystal. The Boson peak is thus 
related to the transverse acoustic singularity of the underlying crystal 
structure. Transverse modes have also been firmly connected to the 
Boson peak in other works 43,56 . Indeed, our studies of transverse and 
longitudinal correlation functions suggest that low-frequency trans- 
verse phonons contribute the most to the Boson peak intensity in the 
range of wave numbers where both the longitudinal and transverse 
phonons show a solid-like response over large length scales, namely 
emergence of a high-frequency peak in both longitudinal and trans- 
verse spectra at cu = 260 cnr 1 for k > 0.5 A -1 . The appearance of this 
high-frequency excitation associated with four-coordinated mole- 
cules in longitudinal and transverse spectra coincides with a pile 
up of intensities in the Boson peak regime. 

To conclude, our results indicate that liquid water displays a Boson 
peak below the Widom line temperature T w . Both the onset temper- 
ature and energy position are similar to what has been observed 
experimentally for confined water 8,9 . We find that as the liquid 
crosses over to a low-density-like liquid structure below T w the 
low-frequency dynamics of the liquid changes to resemble that of 
the underlying crystal, ice Ih. The Boson peak in supercooled water is 
thus a manifestation of the transverse acoustic singularity of the 
crystal and may therefore be a general phenomenon in tetrahedral 
liquids showing a liquid-liquid phase transition. 

Methods 

We simulate TIP4P/2005 36 water for a range of temperatures at atmospheric pressure. 
Equilibration is first performed in the NPT ensemble, using the Nose-Hoover ther- 
mostat and Parrinello-Rahman barostat to attain constant temperature and pressure. 
The equilibrated densities are then used in equilibration NVT runs performed over 
multiple structural relaxation times, after which we switch to the NVE ensemble to 
compute the relevant dynamical quantities. The equations of motion are integrated 
with a time step of 0.2-1.0 fs, depending on the observed energy conservation. Most 
simulations are performed using N = 512 molecules, but to quantify the fmite-size 
effects we simulate larger systems up to N = 45,000. Simulation temperatures 
between 210 and 260 K at P = 1 atm were chosen so that the system crosses T w , the 
temperature where maxima in response functions are observed. To confirm the 
connection between the Widom line T w and the onset temperature of the Boson peak 
we repeat the above simulation protocol for the TIP5P water model 47 at temperatures 
between 240 and 270 K (see Supplementary Information). 



We perform non-equilibrium simulations of TIP4P/2005 at 100 and 150 K by 
rapidly cooling from 210 K with a cooling rate of 2'10 10 K/s, and then switching to 
the NVE ensemble to calculate dynamical properties. 

For the calculation of the participation ratio and the function A max (r), equilibrium 
configurations at a given temperature were quenched to obtain the configurations 
corresponding to the nearest local minimum of the system. Then we calculate 
eigenmodes and eigenvalues corresponding to vibrational modes about this local 
energy minimum by diagonalizing the Hessian matrix with respect to the generalized 
coordinates. We use the flexible version of the model and hence all degrees of freedom 
to calculate the Hessian (see Refs. 54,55,57 for a more detailed explanation of the 
formalism and method). 
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